//#define MULTI_THREAD
#include <process.h>
#include <windows.h>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <vector>

#if __APPLE__
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#include <GL/glext.h>
#endif

#include <Cg/cg.h>
#include <Cg/cgGL.h>

#include "PerlinNoise.h"
#include "TGAImg.h"


#define		TERRAIN_SIZE		256
#define		NOISE_LEVEL			3
#define		NOISE_VOL_W			128			// Width of the noise volume
#define		NOISE_VOL_H			128			// Height of the noise volume
#define		NOISE_VOL_D			128			// Depth of the noise volume
#define		BYTES_PER_TEXEL		1			// 3 unsigned bytes (RGB) for each texel in the volume
#define		PI					3.1415927
#define		TEXTURE_SCALE		1
#define		NUM_OF_TEXTURES		5
#define		NUM_OF_OCTAVES		4
#define		DEFAULT_CAMERA_X	-50.0f
#define		DEFAULT_CAMERA_Y	75.0f
#define		DEFAULT_CAMERA_Z	-150.0f
#define		DISTANCE			sqrt(pow(DEFAULT_CAMERA_X, 2) + pow(DEFAULT_CAMERA_Z, 2))

/* Multithread communication flag */
bool terrainCreated = false;

/* For 3D texture use */
#if __APPLE__
/* Do nothing */
#else
PFNGLTEXIMAGE3DPROC glTexImage3D;
#endif

GLuint noiseTex[NUM_OF_OCTAVES];
unsigned char *Tex3D = NULL;

/* Declare Cg specific variables and functions */
static CGcontext	myCgContext;
static CGprofile	myCgVertexProfile, myCgFragmentProfile;
static CGprogram	myCgVertexProgram, myCgFragmentProgram;
static CGparameter	light, modelView, modelViewMatrix, decal[NUM_OF_TEXTURES], noises[NUM_OF_OCTAVES], fog, animation;
static const char  *myProgramName = "Terrain Vertex",
                   *myVertexProgramFileName = "../Cg/vcolor.cg",
				   *myFragmentProgramFileName = "../Cg/ftexture.cg",
			       *myVertexProgramName = "vmain",
				   *myFragmentProgramName = "fmain";
static void			checkForCgError(const char *);

GLfloat terrain[TERRAIN_SIZE][TERRAIN_SIZE][3];					// Our terrain mesh
GLfloat terrain_norms[TERRAIN_SIZE][TERRAIN_SIZE][3];			// Normal array of our terrain mesh

GLfloat vertices[TERRAIN_SIZE * TERRAIN_SIZE * 3];
GLfloat normals[TERRAIN_SIZE * TERRAIN_SIZE * 3];
GLfloat texCoords[TERRAIN_SIZE * TERRAIN_SIZE * 2];

GLuint	indices[(TERRAIN_SIZE - 1) * (TERRAIN_SIZE - 1) * 2];

GLuint	grassText, sandText, rockText, dirtText, noise1Text;

GLfloat	angle = asin(- DEFAULT_CAMERA_X / DISTANCE) + PI;		// For camera movement
GLfloat zoom = 0;

GLfloat light_position[] = {-1.0, 70.0, -1.0, 0.0};
GLfloat animation_coef[] = {0.12423452, 0.5785673567, 0.4526362, 0.8345215};

GLuint LoadTexture(char *);
GLuint Load3DTexture(int, float);

static int rotating = 0;		// Rotation flag

/* For performance test */
clock_t start, finish;
double duration;

/* For FPS calculation */
clock_t startClock = 0;
int frameCount = 0;

void CrossProduct(float *, const float *, const float *);
void Normalize(float *);
void GenerateTerrainMesh(PVOID param);
void CgSetup();

inline void SetupIndices() {
	int i = 0;
	for (int x = 0; x < TERRAIN_SIZE - 1; x++) {
		for (int z = 0; z < TERRAIN_SIZE - 1; z++) {
			indices[i++] = x * TERRAIN_SIZE + z;
			indices[i++] = (x + 1) * TERRAIN_SIZE + z;
		}
	}
}

bool initialize() {
	
	frameCount++;
	
	/* Start of initialization */
	glClearColor(1.0f, 1.0f, 1.0f, 1.0f);								// White Background
	glClearDepth(1.0f);													// Depth Buffer Setup
	glDepthFunc(GL_LEQUAL);												// The Type Of Depth Testing
	glEnable(GL_DEPTH_TEST);											// Enable Depth Testing
	glEnable(GL_TEXTURE_2D);
	glEnable(GL_TEXTURE_3D);
	glEnable(GL_CULL_FACE);
	glCullFace(GL_FRONT);	
	glShadeModel(GL_SMOOTH);											// Select Smooth Shading
	glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);					// Set Perspective Calculations To Most Accurate
#ifdef MULTI_THREAD
	_beginthread(GenerateTerrainMesh, 0, NULL);
#else
	GenerateTerrainMesh(NULL);
#endif
	SetupIndices();
	init();																// Initialize the generation of Perlin Noise
	return true;
}

void display() {
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	
	/* Camera position */
	gluLookAt((DISTANCE + zoom) * sin(angle), DEFAULT_CAMERA_Y + zoom * DEFAULT_CAMERA_Y / DISTANCE, (DISTANCE + zoom) * cos(angle), 0.0f, 0.0f, 0.0f, 0, 1, 0);
	
	glColor3f(0.5f, 1.0f, 0.5f);
	
	glEnableClientState(GL_VERTEX_ARRAY);
	glVertexPointer(3, GL_FLOAT, 0, vertices);

#if 0	/* Deferred shading */
	glEnable(GL_POLYGON_OFFSET_FILL);
	glPolygonOffset(1, 6);
	for (int i = 0; i < (TERRAIN_SIZE - 1); i++)
		glDrawElements(GL_QUAD_STRIP, (TERRAIN_SIZE - 1) * 2, GL_UNSIGNED_INT, indices + i * (TERRAIN_SIZE - 1) * 2);
	glPolygonOffset(0.0f, 0.0f);
	glDisable(GL_POLYGON_OFFSET_FILL);
#endif
	
	glEnableClientState(GL_NORMAL_ARRAY);
	glEnableClientState(GL_TEXTURE_COORD_ARRAY);
	
	glNormalPointer(GL_FLOAT, 0, normals);
	glTexCoordPointer(2, GL_FLOAT, 0, texCoords);
	
	cgGLEnableProfile(myCgVertexProfile);
	cgGLEnableProfile(myCgFragmentProfile);

	cgGLBindProgram(myCgVertexProgram);
	cgGLBindProgram(myCgFragmentProgram);

	for (int i = 0; i < NUM_OF_TEXTURES; i++)
		cgGLEnableTextureParameter(decal[i]);

	for (int i = 0; i < NUM_OF_OCTAVES; i++)
		cgGLEnableTextureParameter(noises[i]);
		
	cgGLSetParameter4fv(animation, animation_coef);
	animation_coef[0] += 0.00045134531;
	animation_coef[1] += 0.000193124151;
	animation_coef[2] += 0.000653451341;
	animation_coef[3] += 0.000214235521;
	cgGLSetParameter3fv(light, light_position);
	cgGLSetStateMatrixParameter(modelViewMatrix, CG_GL_MODELVIEW_PROJECTION_MATRIX, CG_GL_MATRIX_IDENTITY);
	cgGLSetStateMatrixParameter(modelView, CG_GL_MODELVIEW_MATRIX, CG_GL_MATRIX_IDENTITY);

	start = clock();
#if 0	/* for performance comparison */
	for (int x = 0; x < TERRAIN_SIZE - 1; x++) {
		glBegin(GL_QUAD_STRIP);
		for (int z = 0; z < TERRAIN_SIZE - 1; z++) {
			glNormal3fv(terrain_norms[x][z]);
			glTexCoord2f(((float)x / TERRAIN_SIZE) * TEXTURE_SCALE, ((float)z / TERRAIN_SIZE) * TEXTURE_SCALE);
			glVertex3f(terrain[x][z][0], terrain[x][z][1], terrain[x][z][2]);
			glNormal3fv(terrain_norms[x+1][z]);
			glTexCoord2f(((float)(x + 1) / TERRAIN_SIZE) * TEXTURE_SCALE, ((float)z / TERRAIN_SIZE) * TEXTURE_SCALE);
			glVertex3f(terrain[x+1][z][0], terrain[x+1][z][1], terrain[x+1][z][2]);
		}
		glEnd();
	}	
#else
	for (int i = 0; i < (TERRAIN_SIZE - 1); i++)
		glDrawElements(GL_QUAD_STRIP, (TERRAIN_SIZE - 1) * 2, GL_UNSIGNED_INT, indices + i * (TERRAIN_SIZE - 1) * 2);
#endif
	finish = clock();
	duration = (double)(finish - start) / CLOCKS_PER_SEC;
	//fprintf(stderr, "Time spent drawing the mesh: %.3f second(s)\n", duration);

	cgGLDisableProfile(myCgVertexProfile);
	cgGLDisableProfile(myCgFragmentProfile);
	for (int i = 0; i < NUM_OF_TEXTURES; i++)
		cgGLDisableTextureParameter(decal[i]);
	
	for (int i = 0; i < NUM_OF_OCTAVES; i++)
		cgGLDisableTextureParameter(noises[i]);

	glDisableClientState(GL_NORMAL_ARRAY);
	glDisableClientState(GL_VERTEX_ARRAY);
	glDisableClientState(GL_TEXTURE_COORD_ARRAY);
	
	glutSwapBuffers();

	frameCount++;
}

void reshape(int width, int height) {
	glViewport(0, 0, (GLsizei)(width), (GLsizei)(height));				// Reset The Current Viewport
	glMatrixMode(GL_PROJECTION);										// Select The Projection Matrix
	glLoadIdentity();													// Reset The Projection Matrix
	gluPerspective (45.0f, (GLfloat)(width)/(GLfloat)(height),			// Calculate The Aspect Ratio Of The Window
					0.1f, 1000.0f);		
	glMatrixMode (GL_MODELVIEW);										// Select The Modelview Matrix
	glLoadIdentity();													// Reset The Modelview Matrix
}

void idle(void) {
	
	if (rotating) {
		angle = angle + 0.003;
		if (angle > (2 * PI))
			angle = angle - 2 * PI;
	}
	glutPostRedisplay();
	float seconds = (float)(clock() - startClock) / (float)CLOCKS_PER_SEC;
	float fps = frameCount / seconds;
	printf("\rAverage FPS: %2.2f\t\t\t", fps);
}

void keyboard(unsigned char c, int x, int y) {
	static GLint mode;
	static float fogMode;

	switch (c) {
		case '=':
			zoom -= 10;
			break;
		case '-':
			zoom += 10;
			break;
		case '0':
			zoom = 0;
			break;
		case 'r':
		case 'R':
			rotating = !rotating;
			if (rotating) {
				frameCount = 0;
				startClock = clock();
			}
			break;
		case 'l':
		case 'L':
			glGetIntegerv(GL_POLYGON_MODE, &mode);
			if (mode == GL_FILL)
				glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
			else glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
			break;
		case 'f':
		case 'F':
			cgGLGetParameter1f(fog, &fogMode);
			fogMode++;
			if (fogMode > 2) fogMode = 0;
			cgGLSetParameter1f(fog, fogMode);
			break;
		case ' ':
			//GenerateTerrainMesh();
			break;
		case 27:	/* Esc */
			cgDestroyProgram(myCgVertexProgram);
			cgDestroyProgram(myCgFragmentProgram);
			cgDestroyContext(myCgContext);
			exit(0);
			break;
	}
	glutPostRedisplay();
}

int main(int argc, char** argv) {

	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(1280, 768);
	glutInitWindowPosition(100, 100);
	glutCreateWindow("Fogged Terrain");
	start = clock();
	initialize();
	glutIdleFunc(idle);
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);
#if __APPLE__
/* Do nothing */
#else
	glTexImage3D = (PFNGLTEXIMAGE3DPROC)wglGetProcAddress("glTexImage3D");
	if (glTexImage3D == NULL) {
		fprintf(stderr, "Error in line %d: Couldn't load glTexImage3D function. Aborting.\n", __LINE__);
		return -1;
	}
#endif

	/* Load textures */
	grassText = LoadTexture("../grass.tga");
	sandText = LoadTexture("../sand.tga");
	rockText = LoadTexture("../rock.tga");
	dirtText = LoadTexture("../dirt.tga");
	noise1Text = LoadTexture("../noise1.tga");
	for (int i = 0; i < NUM_OF_OCTAVES; i++)
		noiseTex[i] = Load3DTexture(i, (float)1 / (i + 1));
		
	CgSetup();

#ifdef MULTI_THREAD
	while(!terrainCreated) ;
#endif
	finish = clock();
	fprintf(stderr, "Time used: %f s\n", float(finish - start) / 1000);

	glutMainLoop();
		
	return 0;
}

void GenerateTerrainMesh(PVOID param) {
	std::vector<float *> *normal_buffer = new std::vector<float *>[TERRAIN_SIZE * TERRAIN_SIZE];

	/* Create a flat mesh */
	for (int x = 0; x < TERRAIN_SIZE; x++) {
		for (int z = 0; z < TERRAIN_SIZE; z++) {
			terrain[x][z][0] = (float) (TERRAIN_SIZE / 2) - x;						// Centers the mesh around the origin
			terrain[x][z][1] = 0.0;
			terrain[x][z][2] = (float) (TERRAIN_SIZE / 2) - z;						// Centers the mesh around the origin
		}
	}
	
	/* Determine global shape */
	srand(clock());
	float seed1 = 0.3242352352 + (float)(rand() / 1678);
	float seed2 = 0.5363145 + (float)(rand() / 1042);
	for (int x = 0; x < TERRAIN_SIZE; x++) {
		for (int z = 0; z < TERRAIN_SIZE; z++) {
			terrain[x][z][1] += pnoise(cos(terrain[x][z][2] / TERRAIN_SIZE * PI * 2) + sin(terrain[x][z][0] / TERRAIN_SIZE * PI * 2), 
									   0.3242352352,
									   sin(terrain[x][z][2] / TERRAIN_SIZE * PI * 2) + cos(terrain[x][z][0] / TERRAIN_SIZE * PI * 2))
								* TERRAIN_SIZE / 8;
		}
	}
	
	/* Add some levels of details */
	float noises[NOISE_LEVEL] = {1.131903429420923, 2.39482095842, 3.498058343};
	for(int k = 0; k < NOISE_LEVEL; k++) {
		for (int x = 0; x < TERRAIN_SIZE; x++) {
			for (int z = 0; z < TERRAIN_SIZE; z++) {
				terrain[x][z][1] += pnoise(terrain[x][z][0] / TERRAIN_SIZE * noises[k], terrain[x][z][2] / TERRAIN_SIZE * noises[k], 0.5363145) * TERRAIN_SIZE / (noises[k] * 2);
			}
		}
	}

	/* Calculate the normals */
	for (int x = 0; x < TERRAIN_SIZE - 1; x++)
		for (int z = 0; z < TERRAIN_SIZE - 1; z++) {
			float v1[3], v2[3];
			v1[0] = terrain[x][z][0] - terrain[x+1][z][0];
			v1[1] = terrain[x][z][1] - terrain[x+1][z][1];
			v1[2] = terrain[x][z][2] - terrain[x+1][z][2];
			v2[0] = terrain[x+1][z][0] - terrain[x][z+1][0];
			v2[1] = terrain[x+1][z][1] - terrain[x][z+1][1];
			v2[2] = terrain[x+1][z][2] - terrain[x][z+1][2];
			CrossProduct(terrain_norms[x][z], v2, v1);
			Normalize(terrain_norms[x][z]);
#if 1		/* Together with next one, calculate surface normals */
			normal_buffer[x * TERRAIN_SIZE + z].push_back(terrain_norms[x][z]);
			normal_buffer[(x + 1) * TERRAIN_SIZE + z].push_back(terrain_norms[x][z]);
			normal_buffer[x * TERRAIN_SIZE + z + 1].push_back(terrain_norms[x][z]);
			normal_buffer[(x + 1) * TERRAIN_SIZE + z + 1].push_back(terrain_norms[x][z]);
#endif
		}
#if 1
	for (int x = 0; x < TERRAIN_SIZE - 1; x++)
		for (int z = 0; z < TERRAIN_SIZE - 1; z++) {
			terrain_norms[x][z][0] = terrain_norms[x][z][1] = terrain_norms[x][z][2] = 0;
			for (size_t i = 0; i < normal_buffer[x * TERRAIN_SIZE + z].size(); i++) {
				terrain_norms[x][z][0] += normal_buffer[x * TERRAIN_SIZE + z][i][0];
				terrain_norms[x][z][1] += normal_buffer[x * TERRAIN_SIZE + z][i][1];
				terrain_norms[x][z][2] += normal_buffer[x * TERRAIN_SIZE + z][i][2];
			}
			terrain_norms[x][z][0] /= normal_buffer[x * TERRAIN_SIZE + z].size();
			terrain_norms[x][z][1] /= normal_buffer[x * TERRAIN_SIZE + z].size();
			terrain_norms[x][z][2] /= normal_buffer[x * TERRAIN_SIZE + z].size();
		}
#endif
	
	/* Establish Vertex Arrays */
	int i = 0, j = 0, k = 0, l = 0;
	for (int x = 0; x < TERRAIN_SIZE; x++) {
		for (int z = 0; z < TERRAIN_SIZE; z++) {
 			vertices[i++] = terrain[x][z][0];
  			vertices[i++] = terrain[x][z][1];
			vertices[i++] = terrain[x][z][2];
			normals[j++] = terrain_norms[x][z][0];
			normals[j++] = terrain_norms[x][z][1];
			normals[j++] = terrain_norms[x][z][2];
			texCoords[k++] = ((float)x / TERRAIN_SIZE) * TEXTURE_SCALE;
			texCoords[k++] = ((float)z / TERRAIN_SIZE) * TEXTURE_SCALE;
			//texCoords3D[l++] = ((float)x / TERRAIN_SIZE) * TEXTURE_SCALE;
			//texCoords3D[l++] = ((float)terrain[x][z][1] / TERRAIN_SIZE) * TEXTURE_SCALE;
			//texCoords3D[l++] = ((float)z / TERRAIN_SIZE) * TEXTURE_SCALE;
		}
	}
#ifdef MULTI_THREAD
	terrainCreated = true;
#endif
}

void CrossProduct(float *result, const float *a, const float *b) {
	result[0] = a[1] * b[2] - a[2] * b[1];
	result[1] = a[2] * b[0] - a[0] * b[2];
	result[2] = a[0] * b[1] - a[1] * b[0];
}

void Normalize(float *vector) {
	double length;
	length = sqrt(pow(vector[0],2) + pow(vector[1],2) + pow(vector[2],2));
	vector[0] /= length;
	vector[1] /= length;
	vector[2] /= length;
}

void CgSetup() {
	myCgContext = cgCreateContext();
	checkForCgError("Creating context...");
	
	myCgVertexProfile = cgGLGetLatestProfile(CG_GL_VERTEX);
	cgGLSetOptimalOptions(myCgVertexProfile);
	checkForCgError("Selecting vertex profile...");

	myCgFragmentProfile = cgGLGetLatestProfile(CG_GL_FRAGMENT);
	checkForCgError("Selecting fragment profile...");

	myCgVertexProgram = cgCreateProgramFromFile(myCgContext,              /* Cg runtime context */
												CG_SOURCE,                /* Program in human-readable form */
												myVertexProgramFileName,  /* Name of file containing program */
												myCgVertexProfile,        /* Profile: OpenGL ARB vertex program */
												myVertexProgramName,      /* Entry function name */
												NULL);                    /* No extra compiler options */
	checkForCgError("Creating vertex program from file...");

	myCgFragmentProgram = cgCreateProgramFromFile(myCgContext,                /* Cg runtime context */
												  CG_SOURCE,                  /* Program in human-readable form */
												  myFragmentProgramFileName,  /* Name of file containing program */
												  myCgFragmentProfile,        /* Profile: OpenGL ARB vertex program */
												  myFragmentProgramName,      /* Entry function name */
												  NULL);                      /* No extra compiler options */
	checkForCgError("Creating fragment program from file...");

	cgGLLoadProgram(myCgVertexProgram);
	checkForCgError("Loading vertex program...");

	cgGLLoadProgram(myCgFragmentProgram);
	checkForCgError("Loading fragment program...");
	
	light			= cgGetNamedParameter(myCgVertexProgram, "IN.light");
	modelView		= cgGetNamedParameter(myCgVertexProgram, "ModelView");
	modelViewMatrix	= cgGetNamedParameter(myCgVertexProgram, "ModelViewProj");
	animation		= cgGetNamedParameter(myCgVertexProgram, "Animation");
	fog				= cgGetNamedParameter(myCgFragmentProgram, "fog");
	decal[0]		= cgGetNamedParameter(myCgFragmentProgram, "decal0");
	decal[1]		= cgGetNamedParameter(myCgFragmentProgram, "decal1");
	decal[2]		= cgGetNamedParameter(myCgFragmentProgram, "decal2");
	decal[3]		= cgGetNamedParameter(myCgFragmentProgram, "decal3");
	decal[4]		= cgGetNamedParameter(myCgFragmentProgram, "noise1");
	noises[0]		= cgGetNamedParameter(myCgFragmentProgram, "noiseTex0");
	noises[1]		= cgGetNamedParameter(myCgFragmentProgram, "noiseTex1");
	noises[2]		= cgGetNamedParameter(myCgFragmentProgram, "noiseTex2");
	noises[3]		= cgGetNamedParameter(myCgFragmentProgram, "noiseTex3");

	cgGLSetParameter1f(fog, 0);
	cgGLSetTextureParameter(decal[0], grassText);
	cgGLSetTextureParameter(decal[1], sandText);
	cgGLSetTextureParameter(decal[2], rockText);
	cgGLSetTextureParameter(decal[3], dirtText);
	cgGLSetTextureParameter(decal[4], noise1Text);
	cgGLSetTextureParameter(noises[0], noiseTex[0]);
	cgGLSetTextureParameter(noises[1], noiseTex[1]);
	cgGLSetTextureParameter(noises[2], noiseTex[2]);
	cgGLSetTextureParameter(noises[3], noiseTex[3]);
	checkForCgError("Setting decal 2D & 3D textures...");
}

static void checkForCgError(const char *situation) {
	CGerror error;
	const char *string = cgGetLastErrorString(&error);

	if (error != CG_NO_ERROR) {
		printf("%s: %s: %s\n", myProgramName, situation, string);
		if (error == CG_COMPILER_ERROR) {
			printf("%s\n", cgGetLastListing(myCgContext));
		}
		exit(1);
	}
}

GLuint LoadTexture(char *TexName) {
	TGAImg Img;								// Image loader
	GLuint Texture;

	/* Load our Texture */
	if(Img.Load(TexName) != IMG_OK)
		return -1;

	glGenTextures(1,&Texture);				// Allocate space for texture
	glBindTexture(GL_TEXTURE_2D,Texture);	// Set our Tex handle as current

	/* Create the texture */
	if(Img.GetBPP()==24)
		glTexImage2D(GL_TEXTURE_2D,0,3,Img.GetWidth(),Img.GetHeight(),0,
					  GL_RGB,GL_UNSIGNED_BYTE,Img.GetImg());
	else if(Img.GetBPP()==32)
		glTexImage2D(GL_TEXTURE_2D,0,4,Img.GetWidth(),Img.GetHeight(),0,
					  GL_RGBA,GL_UNSIGNED_BYTE,Img.GetImg());
	else
		return -1;

	/* Specify filtering and edge actions */
	glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_LINEAR);
	glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);
	glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S,GL_REPEAT);
	glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T,GL_REPEAT);

	return Texture;
}

GLuint Load3DTexture(int coef, float scale) {
	GLuint Texture;
	int w, h, d;
	w = h = d = (int)pow(2.0, coef + 3);

	Tex3D = new unsigned char[w * h * d * BYTES_PER_TEXEL];
	if (Tex3D == NULL) {
		fprintf(stderr, "Load3DTexture: Fail to allocate memory for the 3D texture.\n");
		return -1;
	}
	
	for (int r = 0; r < d; r++)
		for (int s = 0; s < w; s++) {
			for (int t = 0; t < h; t++) {
				int index = BYTES_PER_TEXEL * (s * w + t) + w * h * r * BYTES_PER_TEXEL;
				double gray = pnoise(fmod((float)r / d, scale) * 1.14135136, 
									 fmod((float)s / w, scale) * 1.23524624,
									 fmod((float)t / h, scale) * 1.26264266)
								   * 128;
				Tex3D[index] = gray;
				//if (r == 0 && coef == 1)
				//	fprintf(stdout, "%f ", gray);
			}
			//if (r == 0 && coef == 1)
			//	fprintf(stdout, "\n");
		}

	glGenTextures(1, &Texture);
	glBindTexture(GL_TEXTURE_3D, Texture);
	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_REPEAT);
	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_REPEAT);
	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_REPEAT);

	glTexImage3D(GL_TEXTURE_3D, 0, 1, w, h, d, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, Tex3D);

	if (glGetError() != GL_NO_ERROR)
		fprintf(stderr, "OpenGL Error!\n");
	
	delete [] Tex3D;
	
	return Texture;
}